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Abstract 

In the analysis of an interface crack between dissimilar elastic materials, the mode of 
crack extension is typically not unique, due to oscillatory behavior of near-tip stresses and 
displacements. This behavior currently limits the applicability of interfacial fracture mechanics 
as a means to predict delamination in layered materials. The Virtual Crack Closure Technique 
(VCCT) is a method used to extract mode I and mode II energy release rates from numerical 
fracture solutions. The mode of crack extension extracted from an oscillatory solution using the 
VCCT is not unique due to the dependence of mode on the virtual crack extension length, A. 

In this work, a method is presented for using the VCCT to extract A-independent crack 
extension modes for the case of an interface crack between two in-plane orthotropic materials. 
The method does not involve altering the analysis to eliminate its oscillatory behavior. Instead, it 
is argued that physically reasonable, A-independent modes of crack extension can be extracted 
from oscillatory solutions. Knowledge of near-tip Helds is used to determine the explicit A 
dependence of energy release rate parameters. Energy release rates are then defined that are 
separated from the oscillatory dependence on A. A modified VCCT using these energy release 
rate definitions is applied to results from finite element analyses, showing that A-independent 
modes of crack extension result. The modified technique has potential as a consistent method for 
extracting crack extension modes from numerical solutions. The A-independent modes extracted 
using this technique can also serve as guides for testing the convergence of finite element 
models. Direct applications of this work include the analysis of planar composite delamination 
problems, where plies or debonded laminates are modeled as in-plane orthotropic materials. 

Introduction 

The problem of delamination in composite materials Jias been studied extensively in the 
literature. From this work, a methodology for analyzing delamination problems has emerged, 
where the problem is treated within the framework of fracture mechanics. In order to predict 
delamination resistance, the problem of an existing interfacial delamination crack is modeled, 
and an applied stress intensity factor, K, or energy release rate, G, is calculated. The applied K 
or G is then compared to a mode-dependent critical value of K or G from experiment. The work 



of Wang and Crossman (1980), O’Brien (1982) and Wang (1982) was some of the earliest to use 
this type of approach. One advantage of a fracture-based approach is that finite energy release 
rates result, even for limiting cases where the length of the delamination crack approaches zero. 
In contrast, if a fully bonded interface is modeled, singular stresses result at the intersection of 
the interface and the free edge. This makes prediction of delamination through use of a fully 
bonded model problematic. 

The mode dependence of critical K or G values requires that an applied stress intensity 
factor be compared to a measured interfacial toughness for the same mode of crack extension. 
Commonly used interfacial toughness tests include the double cantilever beam test, which yields 
a pure mode I interfacial toughness and the end notched flexure test, which yields a toughness 
close to that for pure mode II (see Fig. 1). Other tests have been proposed to measure toughesses 
under various mixed mode conditions e.g. the work of Reeder and Crews (1990). Composite 
interfacial toughness tests are typically performed on unidirectional 0° laminates. Because the 
elastic properties of adjacent plies are the same, the mode of crack extension associated with 
such tests can be unambiguously determined. In contrast, in most applications, delaminations 
occur between plies of different orientation. In such cases, modeling the delamination as a crack 
along a distinct interface between plies typically results in near-tip stresses and displacements 
that oscillate as the tip is approached. As a result, in the interfacial fracture analysis of most 
composite applications, the mode of crack extension is a function of distance from the crack tip 
and is not uniquely defined. Under such conditions, a direct match of mode between the 
application and coupon tests cannot readily be made. A method is needed for extracting crack 
extension modes from interfacial fracture analyses that exhibit oscillatory behavior. In this 
study, a method is proposed for extracting physically relevant crack extension modes from 
oscillatory analyses where the materials defining the interface are modeled as in-plane 
orthotropic materials. 

Background 

A number of methods have been proposed for extracting fracture modes in interfacial 
fracture analyses that exhibit oscillatory behavior. Two existing approaches involve altering the 
analysis to eliminate the oscillatory nature of the solution. Once this is done, a unique mode of 
crack extension can be extracted. The first of these methods is often termed the resin interlayer 
method and was originally proposed for analyzing isotropic fracture problems by Atkinson 
(1977). In its most commonly used form, it involves inserting a thin homogeneous layer (the 
resin interlayer) between the layers forming the interface and placing the crack within it. 
Because the crack tip is fully embedded in a homogeneous material, non-oscillatory stresses and 
displacements result. 
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The resin interlayer method has some appeal in modeling delamination in resin matrix 
composites because they have a resin-rich region between adjacent plies. This suggests that the 
inclusion of a resin interlayer in the modeling is physically reasonable. The method has 
disadvantages, however. First, the additional effort required to apply the method acts as a 
deterrent to the method’s widespread use. Significantly more effort is required in finite element 
analyses to mesh out the interlayer region. Another disadvantage of the resin interlayer method 
is the addition of the interlayer thickness and properties as variables in the problem. How 
sensitive the results of the analysis are to relatively small changes in the modeling of the 
interlayer is not clear. A final disadvantage of this method involves its true relation to the 
physical delamination problem on an interlayer scale. As Atkinson (1977) himself notes, placing 
the crack fully within the interlayer ignores the fact that a crack under mixed mode conditions in 
a homogeneous material will not grow in a self-similar manner. In the general case of mixed- 
mode loading, a crack within an interlayer will not remain there as it grows, but will branch to an 
interface. On an interlayer scale, modeling a delamination as the extension of a single crack may 
also be non-physical. The work of Chai (1992) suggests that macroscopically mode II dominated 
crack extension in brittle interlayers occurs by the growth of pure mode I microcracks ahead of 
the tip which then link up with the main crack tip. 

A second method of approaching the interface delamination problem involves changing 
“physically insignificant” properties of one or both layers forming the interface to make the 
oscillatory exponent parameter, e, equal zero (definitions for e are given in the next section). 
This method, outlined for orthotropic analyses in the work by Davidson (1993), typically 
involves changing a Poisson’s ratio of one of the layers forming the interface. It is based on an 
analogous method proposed for isotropic bimaterial fracture problems by He and Hutchinson 
(1989) and Suo and Hutchinson (1989). The key to this method is the identification of 
“physically insignificant” properties that can be changed. For the isotropic bimaterial case, it is 
well-established that most interfacial fracture problems are Very weakly dependent on e or the 
related Dundurs parameter p. Thus in most analyses of isotropic bimaterials, it is possible to 
change properties of one or more of the layers forming the interface to make e=0 without 
significantly altering the physics of the analysis. An analogous role for the parameter e for 
orthotropic and anisotropic interface problems has not been established. Although this method 
has potential, standardized application of it requires formulation of consistent criteria for 
identifying which properties can be altered to make e=0 without significantly changing the 
physical aspects of the model. 

A third method for extracting consistent modes of crack extension from oscillatory 
analyses does not involve altering the interfacial fracture analysis. Instead, it is recognized that 
severe mode perturbations due to oscillatory effects are confined to a region very close to the 


3 



crack tip. Outside of this region, the mode of crack extension is only weakly dependent on 
distance from the crack tip. An argument is made that physically reasonable modes of crack 
extension can be extracted from an unaltered interfacial fracture analysis at a prescribed distance 
from the crack tip that is outside of the zone dominated by e effects. The concepts associated 
with this method were originally detailed by Rice (1988) for the case of an interface crack 
between isotropic layers. This method has been applied by numerous researchers to composite 
delamination problems using the Virtual Crack Closure Technique (VCCT), where a virtual 
crack extension length, A, is used that is large compared to the zone dominated by e effects (see 
O’Brien (1982), for example). This method is appealing because it does not require that the 
analysis itself be altered. It results in reasonable and consistent modes of crack extension if a 
consistent value of A is used from problem to problem. The principal drawbacks of this method 
are related to its application. The oscillatory nature of the interfacial fracture solution places a 
limit on how close to the crack tip the VCCT can be applied. However, as with other “local 
fitting” methods, it is important to apply the VCCT in a region close to the crack tip where the 
near-tip (singular) fields dominate. Use of numerical results to identify a value of A satisfying 
both of these requirements can be difficult, particularly if it is coupled with the task of obtaining 
a finite element solution that is sufficiently refined with respect to extracting crack extension 
mode. 

Any of the three methods outlined thus far can, if carefully applied, provide useful mode 
values for interfacial fracture problems. However, each method currently has drawbacks related 
to the ease and/or consistency with which it can be applied by practitioners. In the current study, 
it is argued that practitioners will not want to alter interfacial fracture analyses (e.g. by changing 
properties or by inserting a resin layer). Altering analyses requires increased effort and standard 
methods for altering analyses have not been defined. As in the third method outlined above, 
attention is focused on extracting a physically reasonable mode of crack extension from the 
existing oscillatory analysis. This is done by using asymptotic near-tip stress and displacement 
relations to define a mode that is separated from oscillatory effects. This definition is proposed 
as a consistent measure of crack extension mode for orthotropic interfacial fracture problems. 

Theory 

Isotropic interfacial fracture 

Before the problem of an interfacial crack between orthotropic layers is addressed, the 
problem of an interfacial crack between isotropic layers will be briefly outlined. Although the 
equations associated with the orthotropic problem are more complicated, the conceptual 
problems associated with the isotropic and orthotropic problems are the same. Expressions given 
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for stress intensity factors and near-tip stress fields follow those used by (among others) Rice 
(1988) and Suo and Hutchinson (1990). 

For a crack along the interface between two isotropic layers (see Fig. 2), the singular 
stresses just ahead of the interfacial crack tip (along the positive x axis) take the following form: 


K 


CTyy + 1CTxy V2^ 


( 1 ) 


The near-tip crack face opening and sliding displacements, 82 and 81 (along the negative x axis), 
are given by 

f w ^ 


S 2 + i5j — 


Q +C 2 


2 [(1 + 2ie)cosh(^e) J 

In (1) and (2) the oscillatory exponent, e, is defined as 
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and for isotropic materials 1 and 2 as designated in Fig. 2, P is the second Dundurs parameter 
defined by 


n ^(*2-1) — A* 2 fa-1) (4) 

H 1 (k 2 +1) + h 2 (k 1 + 1)' 

In (4) (ij (j=l,2) is the material shear modulus and Kj = (3-vj)/(l+vj) (j=l,2) for plane stress and 
Kj = 3-4vj (j=l,2) for plane strain. In (2) Cj= (Kj+l)/|ij (j=l,2). In (1) and (2) K=Ki+iK2 is the 
complex stress intensity factor for interfacial fracture problems, where Arabic subscripts are used 
to differentiate Kj and K2 from classically defined stress intensity factors. K takes the 
dimensional form 

K = Kj + iK 2 = f x (applied stress) x (Vhh" l£ ) 


where f is nondimensional and, in general, a complex function of the material properties and the 
specimen geometry. The parameter h is the dominant length scale or characteristic length of the 
near-tip problem. For a long delamination crack between plies of equal thickness the 
characteristic length is the ply thickness. For plies of unequal thickness the characteristic length 
is the smaller of the two ply thicknesses. 

For an interface crack between isotropic layers, the aforementioned oscillatory behavior 
of the near-tip stresses and displacements is due to the x ie term in (1) and (2). As x approaches 
zero the ratio of a xy to cr yy or 81 to 82 changes rapidly. As a result, the mode of crack extension, 
which is related to these ratios, is not unique. One approach (which has not yet been discussed) 
taken to define a mode of crack extension in isotropic interfacial fracture problems involves 
using the phase angle of the complex K as a measure of mode. The mode of crack extension is 
defined by the angle, y, given by 
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( 6 ) 


y/ = tan 


Im(Kh i£ )' 
Re(Kh i£ ) * 


where \j/ is defined to be independent of the characteristic length, h. Inspection of (1) reveals that 
V represents the relative amounts of normal and shear stresses ahead of the crack tip separated 
from the oscillatory quantity (x/h) lE . This method has been suggested and used by (among 
others) Rice (1988), Hutchinson, Mear and Rice (1987) and Charalambides, Lund, Evans and 
McMeeking (1989) as a means for comparing modes of crack extension for isotropic interfacial 
fracture problems. 


Orthotropic Interfacial Fracture 

Consider the interface crack problem illustrated in Figure 2, however materials 1 and 2 
are now orthotropic with the x-y plane as a plane of material symmetry. Suo (1990) has 
formulated equations for the near tip stresses and displacements for such problems. The stress 
field just ahead of the interfacial crack tip (along the positive x axis) takes the form 

(K, +iK 2 )x l£ 

V 2;rx 
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(7) 


The near-tip crack face opening and sliding displacements, 82 and 81 (along the negative x axis), 
are given by 

2H n (K 1+ iK 2 ,)lxp e) 




I — "btOj — . 
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(8) 


In (7) and (8) the subscript H in the complex quantities Oh and 8 h is used to emphasize that the 
normal stress, a yy and crack face opening displacement, 82, are multiplied by Hy quantities 
defined below. In (7) and (8) the oscillatory exponent, e, is given by 


1 , 

£ = — In 

2n 


1 zl 

l + P 


(9) 


which is the same form as that used in the isotropic case. For the orthotropic case, the parameter 
P is given by 


_ [-i/S n S 2 2 + s [2 ] 2 [^S,,S 22 +s, 2 j i 
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( 10 ) 
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where [ ] 2 and [ denote the enclosed quantities evaluated for materials 2 and 1, respectively. 
The quantities Hu and H22 are functions of the elastic properties of the layers forming the 
interface and are defined by 

H n = +[2nA^JS^l (11) 


H 22 = 2nA"^^S^] ] +[2nA"^ A /S^ 
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( 12 ) 


Nondimensional parameters used to define Hu and H22 are 




P = 


2S 12 + S 66 



Finally, for plane stress problems, the Sy in (10), (11) and (12) are the reduced compliances (i 
and j take on values of 1 through 6). For plane strain problems, the Sy are replaced by Sy’ 


defined by 


S i3 S i3 

O » _ O 13 J3 

^33 


(13) 


The problems associated with oscillating behavior for orthotropic interfacial fracture 
problems are the same as those associated with the isotropic case. The oscillatory behavior of 
the near-tip stresses and displacements comes from the x ie terms in (7) and (8). Additionally, 
however, the Hy quantities in the complex quantities cth and Sh complicate the task of defining a 
mode for the orthotropic problem. Suo (1990) suggests using the phase of K as defined for the 
isotropic case in (6) as a consistent measure of mode for orthotropic interfacial fracture 
problems. Unfortunately, due to the Hy quantity multiplying G 22 in (7), this mode definition 
does not have a clear physical significance as is does for the isotropic problem. It also does not 
simplify to the classical definition of mode for cases where e=0 but Hn*H22- In the next 
section, a more physically meaningful definition of mode is introduced. As has been done in 
isotropic interfacial fracture problems, a mode of crack extension is defined that is separated 
from the non-dimensionalized oscillatory behavior of the solution. However, the mode 
definition presented in the current study is related to the relative amounts of energy released by 
an extending crack due to opening vs. that due to sliding. After this mode definition is presented, 
a method is outlined for extracting it from orthotropic fracture analyses exhibiting oscillatory 
behavior using the virtual crack closure technique. 


Mode Definition for Orthotropic Interfacial Fracture Problems 

For the interface crack geometry shown in Fig. 2, the energy release rate for crack 
extension in the x direction can be expressed as 

G total = J 0 A [cr 22 (x)5 2 (A - x) + C7 12 (X)5,(A - x)] dx, (14) 

where the arguments of Irwin (1957) have been used to express the energy per unit width 
released in propagating the crack a distance A along the x axis as the energy of the stresses ahead 
of the tip acting through the displacements behind the tip. Dividing by A and taking the limit as 
A approaches zero gives the energy release rate. In (14) Gtotal is the total energy release rate for 
the extending crack. The first term in (14) is commonly designated as the mode I component of 
Gtotal. corresponding to the energy released by normal stresses acting through crack face opening 


7 



displacements. The second term in (14) is commonly designated as the mode II component of 
Gtotal. corresponding to the energy released by shear stresses acting through crack face sliding 
displacements. The oscillatory nature of orthotropic interfacial fracture solutions causes the 
mode I and mode II components of Gtotal to he functions of the crack extension length A. No 
limiting value of crack extension mode is reached in the limit as A approaches zero. Although 
the individual components of G are a function of A, Gtotal is A independent. 

In this section, the explicit oscillatory nature of the mode I and mode II components of 
Gtotal (designated respectively as Gi and Gn) is determined for orthotropic interfacial fracture 
problems. Some of the procedures used follow those outlined by Raju, Crews and Aminpour 
(1988) for the isotropic case. Modified energy release rate definitions are then presented that are 
separated from oscillatory quantities and are thus independent of the virtual crack extension 
length, A. The mode of crack extension based on these definitions thus corresponds to the ratio 
of energy release due to crack face sliding vs. that due to crack face opening, separated from the 
oscillatory portion of the solution. As (14) indicates, a definition of mode based on energy 
release rates involves products of stresses ahead of the crack tip and displacements of the crack 
faces behind the tip. Mode definitions could also be formulated based individually on the ratio 
of normal stress vs. shear stress ahead of the crack tip or on the ratio of opening vs. sliding 
displacements of the crack faces behind the tip. For a crack in a homogeneous material, all three 
of these mode definitions (based on energy, stress or displacements) are identical. As (7) and (8) 
show, this is not the case for an interface crack, due to the factor (l+2ie) in the denominator of 
(8). This factor causes a small phase shift to exist between the stresses ahead of the tip and the 
displacements behind the tip. 

Although the magnitude of the phase shift between near-tip stresses and displacements is 
small, its existence requires that some choice be made between the three methods for defining a 
mode of crack extension. For the orthotropic problem, the -use of a mode definition based on 
energy release rates is preferred because of the Hjj factors in the complex quantities an and 8 h- 
The Hjj factors, which cannot be separated from the an and 8 h quantities, cancel in the 
calculation of energy release rates. As a result, a mode definition based on energy release rates 

simplifies to the classical definition of mode for cases where e=0 but Hn^H 22 - 

Because the oscillatory behavior of c H and 8^ are known, the initial step in determining 

the oscillatory behavior of Gi and Gn is to express them in terms of <? H and 8 H . First, define two 
real functions d>i and <I>2 in terms of a H and 8 H 

$1 = Re[j Q A ^ H (x) S h (A - x) dxj (15) 

^2 s ^[J^hOO ^h(A-x) dxj, 
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where the overbar designates the complex conjugate of the quantity below it. Using the energy 
release rate definitions given in (14), Gj, Gn and Gtotal can be expressed in terms of <I>i and d>2 
as follows 

(16) 

^ici — G, +G n = lim^^[4>,] . 

Through the definitions for d>i and d>2 the energy release rates have now been expressed in terms 
of Opj and 8^. The oscillatory behavior of d>i and d>2 can be determined by substituting the 


asymptotic stress and displacement equations (7) and (8) into (15) and integrating. The forms of 
d>l and <X >2 to be integrated are 
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Use of the trigonometric substitution x=Asin 2 ((3) in the integral in the expression for <I>i gives an 
expression in terms of the Beta function. Identities between the Beta and Gamma functions 
allow this expression to be simplified so that d>] is given by 


o,= 


H, 


2 cosh Jte 


KK A, 


(18) 


which is non-oscillatory. Use of the same trigonometric substitution in the integral contained in 
G >2 also results in an expression in terms of the Beta function. The expression for d>2 cannot be 
simplified in the way that d>i can. The function d>2 is thus given by 


* 2AH 11 n I 

0 2 = ^ — Re 

7t cosh 7t£ 


jfcos 2 (/J) (cos (P) sin(/J)) 2,f d/J 


(19) 


where the integral is equivalent to B(z,w) with z=l/2+ie and w=3/2+ie. In (19) the crack 
extension length. A, has been normalized by the characteristic length, h. With respect to the 
near-crack-tip fields, h is the dominant length scale determining the size of the region in which 
the near-tip oscillations dominate. For the case of a long delamination crack, the layer thickness 
(or the smallest layer thickness if the layers forming the interface have unequal thicknesses) is 
the characteristic length. 

By determining the oscillatory behavior of d> 2 , the oscillatory behaviors of the energy 
release rate quantities for orthotropic interfacial fracture problems have been isolated. The 
oscillatory behavior of d >2 can be eliminated by defining a d> 2 ’ as 
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( 20 ) 



•J 0 V H (x) ^(A-x) dx 


which is <J>2 with the oscillatory quantity (A/h) 2iE extracted. Now define Gi, G2 (designated with 
Arabic subscripts to differentiate them from Gi and Gn in (16)) and G to tal as 

(2l) 

= G, + G; = UniT^-[<I>,] . 


In these modified energy release rate definitions, Gtotal is still A independent, however, the Gi 
and G2 portions are also A independent. The relations given in (20) and (21) define energy 
release rate quantities that are independent of the crack extension length, A, for orthotropic 
material oscillatory interfacial fracture models. These relations have been derived by first 
isolating the (A/h) 2ie oscillatory behavior of the quantity d>2 and then defining energy release rate 
components in terms of a quantity d>2’ that has this oscillatory quantity extracted from it. The 
result is mode I and mode II energy release rate components that are separated from the 
oscillatory behavior of the solution. 

The definitions for individual mode components derived in this section are analogous to 
the mode definition discussed in the section on isotropic interfacial fracture problems (see the 
discussion related to (6)). In that case, the mode of crack extension was related to the ratio of 
normal stresses to shear stresses ahead of the crack tip, separated from the oscillatory quantity 
(r/h) ie . The same approach is used in deriving (20) and (21), except that mode components are 
given in terms of energy release rates. A phase angle, \|/, analogous to that defined in (6) can be 
used to designate mode mix. 


y/ = tan 



( 22 ) 


where Gi and G2 are given by (21). 

The basis for the mode definitions derived in this section is that the near-tip oscillatory 
behavior of interfacial fracture solutions is non-physical. It is therefore argued that a mode of 
crack extension defined to be separated from oscillatory quantities in the solution can serve as a 
reasonable and consistent definition for comparing modes between problems. The normalization 
of A by h used in (20) is chosen because h is the length scale determining the size of the region in 
which the near-tip oscillations dominate. It is important to note, however, that this normalization 
is not the only one that could be reasonably argued. For example, one might argue that the zone 
dominated by the near-crack-tip stresses and displacements (often referred to as the K field) will 
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extend to a distance that is a fraction of h. Therefore, a normalization of A by some fraction of h 
in (20) could be proposed. The validity of any proposed normalization (including the one 
proposed in this paper) must eventually be proven by its ability to consistently relate interfacial 
fracture resistances. It is important to recognize, however, that the numerical difference in the 
mode components that would result from any physically reasonable change in normalization 
would be small. Because the magnitude of e is very small for most physical problems (see Suo 
(1990)), the mode change from oscillations over a range of distances from the crack tip of h to 
h/20 is very small. From a practical standpoint, determination of which normalization is the 
“correct” one would be difficult and would likely have minimal bearing on the physical task of 
consistently comparing modes of crack extension between fracture problems. 


Virtual Crack Closure Technique 

The virtual crack closure technique (VCCT) is a method used to extract mode I and mode 
II energy release rate components from numerical (typically finite element) analyses. The 
method, which is detailed by Rybicki and Kanninen (1977), is based on the Irwin representation 
of energy release rate for self-similar crack extension given in (14). The method involves 
expressing the mode I and mode II portions of Gtotal (see (14)) in terms of nodal forces ahead of 
the crack tip and nodal displacements behind the tip. This can be represented by the expressions 


G,=-EfJ 2 G n = — Z F 12 5 1 , 

2 A i 2 i 2 J 2A J nodes i ^ 


(23) 


where F 22 j and Fj 2 j are, respectively, the normal and shear nodal forces per unit thickness over 
distance A ahead of crack tip and § 2 j and Sjj are "corresponding" nodal crack opening and 
sliding displacements over distance A behind the crack tip. Corresponding nodes are designated 
in Figure 3. In the VCCT, A is referred to as the virtual crack extension length. The VCCT is 
applied by using (23) to calculate energy release rates for smaller and smaller values of A, until 
no significant changes in G values occur. In applying the VCCT to oscillatory interface crack 
solutions, only Gtotal reaches a limiting value as A is decreased. The relations in (23) assume that 
nodes are equally spaced ahead of and behind the crack tip, as would be the case if non-singular 
elements are used in the analysis. The technique can also be used with singular crack tip 
elements using relations given in Raju, Crews and Aminpour (1988). 


Modified Virtual Crack Closure Technique 

In this section, a modified VCCT is formulated to allow extraction of the energy release 
rates defined in (21) from a numerical model that exhibits oscillatory behavior. To do this, the 
energy release rate quantities in (21) are equivalently expressed in terms of nodal forces and 
displacements. Define 
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( 24 ) 



where F 22 j and F^j and § 2 j and 5 jj are defined with respect to Fig. 3 as they are for the 
traditional VCCT. With these complex quantities defined, d>i and d> 2 ’ can be expressed as 



and the energy release rate (Gi, G 2 and G to tal) definitions in terms of d>i and $>2 are simply 
those given in (21). 

In applying this modified virtual crack closure technique, the same quantities (nodal 
forces ahead of tip and nodal displacements behind the tip) are used as in the traditional VCCT. 
The quantities are simply re-packaged in the complex forms given in (24). The only additional 
effort required is that needed to calculate the quantities Hu, H 22 and e (defined in (9) through 
(13)) and to perform the complex multiplication indicated in (25). Use of the modified VCCT 
outlined in this section is similar to the existing approach of using the traditional VCCT with a 
virtual crack extension length. A, that is large compared to the zone dominated by e effects (this 
approach was described in the final portion of the Introduction). The goal of using the traditional 
VCCT with a large value of A is to avoid oscillatory effects by simply avoiding the near-tip 
region. The advantage of this modified VCCT is that small values of A can be used while still 
avoiding oscillatory effects because the oscillatory effects are extracted analytically. This allows 
limiting values of energy release rates to be obtained as A approaches zero. 

A relationship exists between the modes of crack extension obtained using the traditional 
and modified VCCT. The mode of crack extension obtained using the modified technique 
corresponds to the mode one would obtain using A=h in the traditional VCCT, if the zone 
dominated by the near-crack-tip fields extended to a distance h from the crack tip. This 
relationship is due to the normalization with respect to h used in the modified VCCT (see (19) 
and (25)). Because a second interface exists at a distance h from the crack tip, in practical 
problems the near-tip zone will not extend to a distance A=h. As a result, direct comparison of 
the two methods based on their predictions for A=h is generally not possible. 

Application 

In the previous section, a modified VCCT has been formulated to extract the energy 
release rate quantities defined in (21). In theory, application of this modified VCCT to numerical 
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solutions will result in A-independent modes of crack extension. This section demonstrates how 
the methods outlined in this study perform in practice. Results of application of the modified 
VCCT to finite element models of two orthotropic interfacial fracture problems are presented. 
An example of the finite element mesh used in this study is shown in Figure 4. The model is 
constructed out of eight-noded plane strain quadrilateral interpolation elements using the finite 
element package ABAQUS. Refined singular and nonsingular element meshes have been used 
near the crack tip. For the nonsingular near-tip mesh, nodes are located at equally spaced r 
intervals from the crack tip. For the singular near-tip mesh, the first elements surrounding the tip 
are quarter-point elements to capture the 1/Vr near-tip strain dependence. Subsequent elements 
are biased toward the crack tip to give higher element resolution there. The density of the near- 
tip mesh was varied to check for convergence; however, for the results presented here, the 
singular and non-singular near- tip meshes consist of 18 rings of elements meshed over a length 
equal to h/2, where h is the smaller thickness of the two plies forming the interface. 

Two problems were analyzed to evaluate the methods proposed in the previous sections. 
The problems analyzed are shown in Figure 5 and are designated as test case #1 and test case #2. 
The two test cases consist of plane strain drop-ply configurations of 0° and 90° 0.005 inch thick 
graphite-epoxy plies. Each test problem has displacements constrained to equal zero on the left 
edge (built-in conditions) and a shear load distributed along the right edge. Test case #1 has been 
previously modeled in work by Wang, et al (1993) on skin-stiffener debond modeling using plate 
elements. Test case #2 was chosen to allow demonstration of the behavior of the modified 
VCCT for both positive and negative values of e. Unidirectional ply properties used in the 
analysis are identical to those used by Wang, et al (1993) and are given in Table 1. Although the 
magnitude of e is the same for both cases, the physical behavior of the problems is very different. 
Test case #2 exhibits more compliant response compared to test case #1. The energy release 
rates for test case #2 are much higher than those for test case #1 and the mode of crack extension 
for the two problems is very different. 

Figures 6 and 7 give plots of the ratios Gn/Gi and G 2 /G 1 as a function of A/h for test case 
#1 and test case #2, respectively. In both figures, the ratio is plotted as calculated using the 
modified and traditional virtual crack closure techniques. In Figure 7, the two points to the far 
left represent values calculated using each virtual crack closure technique over the first element 
of a mesh containing singular crack tip elements. The four pairs of points to the right represent 
values calculated using one, two, three, and four elements from the crack tip of a mesh 
containing non-singular crack-tip elements. The plot in Figure 7 contains only points from a 
model with non-singular crack tip elements. 

Figures 6 and 7 show that application of a traditional VCCT to these two cases using 
small values of A/h results in significant changes in Gn/Gi with A. The positive value of e in test 
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case #1 causes Gn/Gi to decrease with decreases in A in Fig. 6. In Fig. 7, the negative value of e 
in test case #2 causes Gn/Gi to increase with decreases in A. As Figures 6 and 7 indicate, use of 
the modified VCCT results in G 2 /G 1 ratios that are essentially independent of A. Thus the 
modified VCCT appears to work well in allowing extraction of a consistent mode of crack 
extension from an oscillatory solution. Slight variations in the Q/Gi values obtained using only 
the first elements from the crack tip can be seen. These can be attributed to the fact that rapid 
changes in displacements as the crack tip is approached can cause the first element from the tip 
to give less accurate nodal force and displacement results than subsequent elements away from it. 
As mentioned in the previous section, the mode of crack extension extracted using the modified 
VCCT corresponds to that from a traditional VCCT using A=h. The results shown in Figs. 6 and 
7, which show increasing agreement between the two methods as A is increased, qualitatively 
agree with this. Curves for the two methods would continue to approach each other if A were 
increased further, as long as the length A remained within the region dominated by the near-crack 
tip fields. Again, it should be emphasized actual application of the traditional VCCT for values 
of A approaching h would generally give incorrect results because the near-tip fields no longer 
dominate at a distance h from the crack tip. 

The results in Figs. 6 and 7 point out one final advantage of the modified VCCT, related 
to the task of verifying the convergence of a finite element solution. By looking at the modified 
VCCT results in Figs. 6 and 7, one can be reasonably well-assured that the near- tip mesh is 
sufficiently refined. Because the ratio of G 2 /G 1 is constant for the modified VCCT applied over 
a number of elements, it is clear that multiple elements are within the region dominated by the 
near-crack-tip fields. Also, the near-tip elements are correctly picking up the variations in mode 
due to the oscillatory nature of the solution (because when the modified VCCT factors out 
oscillatory quantities a constant value results). Thus the modified VCCT can aid numerical 
modelers in determining if they have a sufficiently refined near-tip mesh, substantially reducing 
the need to compare results from meshes with many levels of refinement to explicitly verify 
convergence. 

Conclusions 

In this study, the problem of extracting consistent modes of crack extension from 
oscillatory near-crack-tip solutions has been addressed. In particular, a mode definition for the 
problem of a crack along an interface between two in-plane orthotropic layers has been 
developed and applied within a modified version of the virtual crack closure technique. 
Knowledge of the near-tip oscillatory stress and displacement fields for an interface crack 
between two in-plane orthotropic materials has been used to determine the explicit oscillatory 
nature of mode I and mode II energy release rate components. Energy release rates have 


14 



subsequently been defined that are separated from oscillatory quantities in the solution. These 
energy release rate definitions have been used to formulate a modified Virtual Crack Closure 
Technique (VCCT) that allows extraction of crack extension modes that are A independent. The 
modified VCCT uses the same quantities (near crack tip nodal displacements and forces) 
extracted from a finite element analysis using the traditional VCCT. The modified technique has 
been applied successfully to two numerical test cases. 

The methods developed in this study can be applied to any in-plane orthotropic 
debonding problem, however the most obvious application is to composite delamination 
problems. The current work is directly applicable to planar composite delamination problems 
where the plies or sublaminates are modeled as in-plane orthotropic layers. The definition for 
mode of crack extension defined in this study has potential as a consistent mode definition to be 
used in such problems. Additionally, the ability to extract A-independent mode values using a 
virtual crack closure technique can aid numerical modelers in determining solution convergence. 

Acknowledgment 

The author gratefully acknowledges the support of the NASA Langley Mechanics of 
Materials Branch and the American Society for Engineering Education through their joint 
Summer Faculty Fellowship Program. The author would also like to thank Kevin O’Brien for 
his thoughts and insight related to this research. 

References 

Atkinson, C. (1977) On stress singularities and interfaces in linear elastic fracture mechanics. 
International Journal of Fracture, 13, 807-820 

Chai, H. (1992) Experimental evaluation of mixed-mode fracture in adhesive bonds. 
Experimental Mechanics, 32, 296-303. 

Charalambides, P.G., Lund, J., Evans, A.G. and McMeeking, R.M. (1989) A test specimen for 
determining the fracture resistance of bimaterial interfaces. Journal of Applied Mechanics, 56, 
77-82. 

Davidson, B.D. (1993) Prediction of energy release rate for edge delamination using a crack tip 
element. Proc. Fifth ASTM Symposium on Composite Materials: Fatigue and Fracture, 
Atlanta, May 4-6. 

He, M.Y., and Hutchinson, J.W. (1989) Kinking of a crack out of an interface. Journal of 
Applied Mechanics, 56, 270-278. 

Hutchinson, J.W., Mear, M.E. and Rice, J.R. (1987) Crack paralleling an interface between 
dissimilar materials. Journal of Applied Mechanics, 53, 828-832. 


15 



Irwin, G.R. (1957) Analysis of stresses and strains near th6 end of a crack traversing a plate. 
Journal of Applied Mechanics, 24, 361-364. 

O'Brien, T.K. (1982) Characterization of delamination onset and growth in a composite 
laminate. Damage in Composite Materials, ASTM STP 775, 140-147. 

Raju, I.S., Crews, J.H., Jr. and Aminpour, M.A. (1988) Convergence of strain energy release 
rate components for edge-delaminated composite laminates. Engineering Fracture 
Mechanics, 30, 383-396. 

Reeder, J.R. and Crews, J.H., Jr. (1990) Mixed-mode bending method for delamination testing. 
AIAA Journal, 28, 1270-1276. 

Rice, J.R. (1988) Elastic fracture mechanics concepts for interfacial cracks. Journal of Applied 
Mechanics, 55, 98-103. 

Rybicki E.F. and Kanninen M. F. (1977) A finite element calculation of stress intensity factors 
by modified crack closure integral. Engineering Fracture Mechanics, 9, 931-938. 

Suo, Z. (1990) Singularities, interfaces and cracks in dissimilar anisotropic media. Proceedings 
of the Royal Society of London, A 427, 331-358. 

Suo, Z. and Hutchinson, J.W. (1989) Sandwich test specimens for measuring interface crack 
toughness. Materials Science and Engineering, A107, 135-143. 

Suo, Z. and Hutchinson, J.W. (1990) Interface crack between two elastic layers. International 
Journal of Fracture, 43, 1-18. 

Wang, A.S.D. and Crossman, F.W. (1980) Initiation and growth of transverse cracks and edge 
delamination in composite laminates parts 1 and 2. Journal of Composite Materials, 
supplemental volume, 71-106. 

Wang, J.T., Raju, I.S., Davila, C.G. and Sleight, D.W. (1993) Computation of strain energy 
release rates for skin-stiffener debonds modeled with plate elements. Proc. 34th 
AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Part 
3, Paper AIAA 93-1501-CP, 1680-1692. 

Wang, S.S. (1982) Fracture mechanics for delamination problems in composite materials. 
Progress in Science and Engineering of Composites, T. Hayashi, K. Kawata, and S. 
Umekawa, Ed., Proc. ICCM-IV, Tokyo, 287-296. 


16 



Table 1 

Graphite-Epoxy Ply Properties Used in the Numerical Test Cases 
(All moduli are in psi) 


E n = 19.5 x 106 E 22 = e 33 = i . 48 x 106 

M -12 = Hl 3 = 0.80 x 10 6 |X 23 = 0.497 x 10 6 

Vi2 = Vi3 = 0.30 V23 = 0.49 


Note: 1 designates the direction along the fiber, 2 designates the in-plane direction normal to 
the fiber and 3 designates the out-of-plane direction normal to the fiber. 

These properties are the same as those used by Wang, et al (1993). 
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Figure 1 Double Cantilever Beam and Edge Notched Flexure Fracture Toughness Specimens 
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Figure 2 Interface Crack Between Two Elastic Layers 
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Figure 3 Virtual Crack Closure Technique Parameters 
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Figure 4 Finite Element Mesh with Non-Singular Crack Tip Elements 
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Figure 6 Energy Release Rate Ratios vs. Normalized Virtual Crack Extension Length, A, for 
Test Case #1 
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Figure 7 Energy Release Rate Ratios vs. Normalized Virtual Crack Extension Length, A, for 
Test Case #2 
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